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ABSTRACT 

In this paper we present a series of models for the deep water cycle on super-Earths experiencing 
plate tectonics. The deep water cycle can be modeled through parameterized convection models 
coupled with a volatile recycling model. The convection of the silicate mantle is linked to the volatile 
cycle through the water-dependent viscosity. Important differences in surface water content are found 
for different parameterizations of convection. Surface oceans are smaller and more persistent for single 
layer convection, rather than convection by boundary layer instability. Smaller planets have initially 
larger oceans but also return that water to the mantle more rapidly than larger planets. Super- 
Earths may therefore be less habitable in their early years than smaller planets, but their habitability 
(assuming stable surface conditions), will persist much longer. 


1. INTRODUCTION 

In the post-Kepler era, the number of known planets 
of Earth and super-Earth size continues to grow. The 
TESS mission, which is a planned survey of the bright¬ 
est and closest stars, is expected t o detect thousands of 
small planets (iRicker et al. l[2015lh The TESS planets, 
by virtue of their closeness and the brightness of their 
stars, will be much easier to follow up with ground-based 
photometric and spectroscopic techniques than the plan¬ 
ets of similar size detected by Kepler. Currently, the 
mini-Neptune GJ 1214 is the only planet in the super- 
Earth mass range (1 — lOAf®) for which atmosp heric 
measurements have been made (e.g..lWilson et al. Il20l4 

iFraine et al. 11201,1 : iBerta et al. It2012f >. With future tele¬ 

scope resources such as GMT and JWST, we will be able 
to characterize even more super-Earths, some of which 
may be in the habitable zones (HZ) of their stars. The 
classical habitable zone is defined as the orbital region in 
which an Earth-like planet can maintain liquid water on 
its surface given a variety of atmospheric compositions. 
However, the Earth’s atmosphere, although it permits 
surface water to remain liquid, does not control the sup¬ 
ply of water on the surface. For that, we must look to the 
mantle. In extending this to super-Earths, it is presently 
unclear what the effect of mass may be in controlling the 
supply of water at the surface. 

Water on the Earth is found in two primary reser¬ 
voirs: (1) the surface oceans and (2) the silicate man¬ 
tle. Water, either in molecular form or as an H + or 
OH - group, can dissolve in silicate minerals. This is 
true even of the so-called nominally anhydrous minerals 
(NAMs), which do not have hydrogen in their chemical 
formula. Furthermore, the solubili ty o f water in sili cate 
minerals increases with p ressure (IKohlstedt et al. 1119961 
iHirschmann et aI7~l 1200511 . 

Geochemical studies of the materials that make up 
Earth’s mantle have long suggested that it may hold any¬ 
where from 1 to 10 times as muc h water as is found in th e 
surface oceans (e.g., Table 1 of iBounama et alHl2001h . 
although more recent work seems to discount larger val- 
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ues in favor of values rangin g from 0. 5 to 2.5 ocea n 
masses (OMs) of water (IHirschmann fe Kohlstedt I2012T) . 
The presence of water in these minerals strongly influ¬ 
ences their material properties, such as melting temper¬ 
atures, rheology, ph ase chan ges, electrical conductivity, 
etc, (see e.g., Hirt.h fc Kohlstedt, 1 119961 : iKarato 1 119901 : 
iLitasov fe Ohtani 1120071 b Ihudies of the effect of water 
on silicate rheology show that the presence of even minor 
amounts of water c an redu ce t he viscosi ty by several or¬ 
ders of magnitude (jKarato fc Wu I il 993t h This can have 
a significant effect on the material flow of the mantle. 

In fact, it has been shown that the abundance of wa¬ 
ter on the Earth’s surface is controlled by the deep wa¬ 
ter/silicate cycle, which is tied to plate tectonics. The 
abundance of water at the surface is a balance between 
the rate of return via volcanic outgassing from the mantle 
at mid-ocean ridges (MORs) and the rate of loss to the 
mantle through so-called ingassing, which is the return of 
water into the deep mantle by the sinking, or subducting, 
of water-rich oceanic seafloor. Much of this water is re¬ 
leased immediately back to the surface through shallow, 
water-induced volcanism. However, a small but signif¬ 
icant fraction of the water contained in the subducting 
oceanic slabs can be transported to deeper levels of the 
mantle. For a detailed review of the exchange rnecha- 
nisms between the surface and mantle, see IHirschmann 1 

(l2006h . 


The deep water cycle on Earth has been stud¬ 
ied through the use of parameterized convection mod¬ 
els incorporating a water-d e pendent viscosity (e.g . 
Mc Govern fe Schubert 119891: IBounama et al. I 1200lc 
Sandu et al. 1 120111: iCrowlev et al. 1 1201 111 . Parameter¬ 
ized convection models are simplified ID models using 
a parameterization derived from more complicated and 
expensive 2D and 3D models of convection for differ¬ 
ent systems (e.g. spherical vs. plane-parallel, heated 
from within vs. below, etc.). The parameterized mod¬ 
els rely primarily on two dimensionless parameters: the 
Rayleigh number, which describes whether a system is 
unstable to convection, and the Nusselt number, which 
compares the convective and conductive heat flows. The 
Rayleigh number depends on the viscosity of the system, 
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which itself is dependent on temperature, pressure and 
water fugacity. The abundance of water in the mantle, 
which helps determine the viscosity, evolves along with 
the mantle temperature. 

Here we will use a parameterized convection model to 
study the deep water cycles of super-Earths. In the 
current era of exoplanet studies, we are still search¬ 
ing for Earth-like planets in Earth-like orbits around 
Sun-like stars, but what we have found and what we 
can soon characterize, are super-Earths (~ 1 — 2i? 0 ,~ 
1 —IOMq). Although super-Earths may have cool enough 
surfaces to retain liquid water, tectonics and mantle 
convection plays an important role in controlling its 
abundance at the surface. The question of whether 
these planets will have plate tectonics has been dis¬ 
cusse d previously in the lit e rature ( e.g. iValencia et, alH 


2007: O’Neill & 

^enardic 2007^ Korenaea 1 

Noack & Breuer 

2013|; Stamenkovic & Breuer 


although the issue has not been settled. Here, we as¬ 
sume plate tectonics as a starting point in order to apply 
a similar water cycle model. The question this work then 
addresses is whether the different pressure regimes inside 
super-Earths affect the deep water cycle, and whether 
surface oceans are tectonically sustainable on these plan¬ 
ets. This has important implications for the potential 
habitability of super-Earths. 

In this paper, we address these questions using a pa¬ 
rameterized thermal evolution model coupled with a wa¬ 
ter cycle model. We give a full description of the model 
in Section 2. In Section 3, we describe results for mod¬ 
els which have either single layer convection or bound¬ 
ary layer convection for a number of super-Earth mod¬ 
els. We also describe the response of the models to rea¬ 
sonable variations in the parameter values. Section 4 
discusses implications for the persistence of oceans on 
super-Earths. 


2. METHODS 

Mantle convection in the terrestrial planets can be 
controlled either by the stability of the who le man¬ 
tle lay er (single layer convection, e.g. iSchubert et al. I 
HOOII) or by the stability of two boundary layers: a 
cold boundary layer at the surface and a hot bound¬ 
ary layer at the inter face of the silic ate ma ntle w ith the 
metallic core (see e.g. iTurcotte fc Schubert If2002l ). Heat 
from either secular cooling or decay of radioactive ma¬ 
terials (or both) is transferred by conduction through 
the boundary layers. The interior region is convective 
and thus approximately adiabatic. See Figure 1 for a 
schematic representation of the thermal profile. Single 
layer convection models typically neglect the heat flux 
from the core, which is small, for simplicity and are 
entirely heated from within by radioactive decay. Nu¬ 
merical simulations are used to determine scaling laws 
relating the heat flux and mean temperature to the de¬ 
gree of convection. These simulations have been done 
for a variety of different boundary conditions (free slip, 
no slip, etc.), and for different material properties (iso- 
visco us, highly temperatur e -dependent viscosity, etc.) 
fe.g.. iHonda fc Iwase 119961 : iDeschamns fc Sotin 1 1200ft 
iDumoulin et al. 112005 ). Choosing the proper scaling re¬ 
lations is therefore important for a successful model. 

The parameter which is of the most impor¬ 
tance in determining the convective behavior is typi¬ 



Figure 1. Schematic of temperature-depth profile. See text for 
details. 


cally the characteristic viscosity. Discussion in the lit¬ 
erature is often contradictory on where to define the 
characteristic viscosity: as a n average value (s e e e.g . 
McGovern fc: Schubert I 119891 iTaiika fc Matsui I Il992t 


Sandu et al. 1 1201 lib or at the base o f a boundary layer 


see e.g. DeschamBS_fe_Sptin_ 120001: IDumoulin et al. 


120051 : IStamenkovic et al. 1 120 IS ). Here we will try both 
types of models, starting with the single-layer models 
typically used with volatile evolution models, which use 
average viscosities. We will then look at a boundary layer 
model. We explore both models here to see the effect on 
the behavior of water, and note that we are not attempt¬ 
ing to reproduce the Earth, but explore two valid, but 
different, models. 

The water cycle model is coupled to the thermal model 
through the viscosity, which is dependent on the w ater 
abundance in the mantle ( see e.g. lMcGovern fc Schubert! 
119891 : ISandu et al. 1 1201 ill . Models of the Earth have 
shown that the vi scosity and mantle temperature create 
a fee dback loop (jSchubert et al. 1 120011: ICrowlev et al. I 
1201 If) . When the mantle is warm, convection is vigor¬ 
ous and the mantle cools quickly. As the temperature 
drops, the viscosity increases, causing the convection to 
become sluggish. Sluggish convection means that less 
heat is removed from the mantle, causing it to heat up. 
The viscosity increases, and the cycle repeats. This cycle 
has been shown to be enhanced by the presence of wa ¬ 
ter (iMcGovern &; Schubertll 1 9891 : [Crowley et al. 11 2 OTTI) . 
ICrowlev et al. (iMlil describe how the water and tem¬ 
perature feedbacks inte ract for temperature an d water- 
dependent viscosities. I Cowan fc Abbotl (1201411 studied 
the effect of sea floor pressure on a steady state model 
of the deep water cycle, without considering the effect of 
planetary thermal evolution. 

In the following sections, we first describe selection of 
the super-Earth model parameters. We then describe 
the parameterized thermal model, followed by the water 
cycle model. Section 3 will discuss the results for these 
models. 


2 . 1 . Super-Earth Models 












































































Oceans on Super-Earths 


3 


Table 1 

Super-Earth Model Parameters. 


Mass 

(M®) 

R p 

(R®) 

Rcore 

( R ®) 

( Pm) 

(kg m“ 3 ) 

(9) 

(m s- 2 ) 

1 

1 

0.547 

4480 

9.8 

2 

1.21 

0.649 

4960 

13.4 

3 

1.34 

0.717 

5470 

16.4 

5 

1.54 

0.814 

5970 

20.7 


Note. — Core mass fraction is held fixed at 0.3259, and 
scaling relations from Valencia et al. (2006) are used to deter¬ 
mine R p and R C ore • See text for details. 


We use the scaling relations of [Valencia et al. 1 (|2006T ) 
to calculate the planetary parameters for planets of 1 , 2 , 
3, and 5 Earth masses (1M© = 5.97 x 10 2 4 kg) . Larger 
planet s are not considered here because lNoack fc Breuer I 
(I2013T) find that the peak likelihood for plate tectonics 
oc curs for planets betwe en 1 — 5M®. The scaling laws 
of [Valencia et al. 1 (l2006h assume a constant core mass 
fraction of 0.3259. The planetary and core radii then 
scale by Ri ~ Ri^(M p /M^) ai , where i = c (core) or 
p (planet), a c = 0.247, a p = 0.27, and values for the 
Earth are indicated by ©. The average mantle density 
(pm) is calculated from the mantle mass and volume. 
The average gravitational acceleration ( g ) is found from 
GM p /R p . Values for these parameters are given in Table 
1. We take a constant water mantle mass fraction of 
1.4 x 10 -3 for the nominal models. This is equivalent 
to 4 ocean-masses (OM) of water for the Earth, where 1 
OM is equal to 1.39 x 10 21 kg H 2 O. We will later explore 
the effect of variable water on the results. 

2 .2. Thermal Evolution Model 

Fo llowing models of Earth’s deep water cycle 

(e.g.lMcGovern fc Schubert~lll989t ISchubert et al. 1120011 : 

[Sandu et al. 11201 lib we will first consider models heated 
only from within (i.e., heat flux q c = 0). These models 
are considered single-layer convection because the whole 
mantle convects. There is a conductive upper thermal 
boundary layer that governs heat loss from the surface. 
For most of the lifetime of the Earth, such models have 
been shown to give g ood fits to the observe d mantle vis¬ 
cosity and heat flux (ISchubert et al. Il2001lh The ther¬ 
mal evolution model requires solution of the mantle heat 
transfer equation: 

PmCpVm ^ ^ = —A s q s + A c q c + V m Q(t) (1) 

where p m is the mantle density, C p is the mantle heat 
capacity, V m is the mantle volume, A Sl A c are the sur¬ 
face area’s of the planet and core, q s and q c are the con¬ 
ducted heat fluxes through the surface and core-mantle 
boundary (CMB), respectively, and Q[t) is the heat pro¬ 
duced by radioactive decay within the mantle. In this 
model, there is no heat flux from the core, so the sec¬ 
ond term on the right side vanishes. The temperature 
modeled in any parameterized convection model is some- 
w hat am bi guous. Here, we will follow the convention of 
iMcGovern fe Schubert I (|l989f ) and take the temperature 
to be the spherically-averaged mantle temperature. We 
can relate this averaged temperature to the temperature 
of the mantle adiabat extrapolated to the surface (i.e. 


the mantle potential temperature T p ) through the equa¬ 
tion: 

(T m ) = e m T p = R3 S _ R3 j* P r 2 T(r)dr (2) 

See lTaiika fc MatsuT ( 199 2) for a derivation of this equa¬ 
tion. We approximate the adiabat as: 

T{r) = T p + T p — Ar (3) 

Cp 

Values of e m for the different planet masses are given in 
Table 5. 

In the mantle, heat is generated by decay of radioactive 
elements. The heat flux from radioactive decay is domi¬ 
nated by 238 U, 235 U, 232 Th, and 40 K. The heat produced 
is calculated from the equation: 

Q(t ) = p m E CiHiexp[\i(4.6 x 10 9 - t)] (4) 


where C) is the mantle concentration of the element by 
mass, Hi is the heat production per unit mass, and A i is 
the decay constant. The decay constants, heat produc¬ 
tion rates, and ab und ances r elativ e to to tal uranium are 
taken from iTurcotte fe Schubert 1 (|2002f ). In this paper, 
we assume that all super-Earths have the same ratios of 
radioactive elements as the present day Earth. The nom¬ 
inal bulk silic ate Earth (i.e., primitive mantle) contains 
~ 21 ppb U ([McDonough fe Sun~lll995lh 
Using boundary layer theory, the heat flux out of the 
mantle is given by: 


Qm — 


(T u - T a ) 


(5) 


where k is the mantle conductivity, S u is the boundary 
layer thickness, and the T„ is the temperature at the base 
of the boundary layer (see Figure 1). T u is calculated 
from (T m ) using the adiabatic temperature profile of the 
mantle. 

The thickness of the upper boundary layer is given by 
the global Rayleigh number: 



/ kt]{T, P)Ra crit \ 
V gap m AT ) 


( 6 ) 


where Z is the thickness of the mantle, Ra is the 
Rayleigh number of the whole mantle, Ra cr n is the crit¬ 
ical Rayleigh number for convective instability, n is the 
mantle thermal diffusivity, r](T, P) is the characteristic 
mantle viscosity, and a is the thermal expansivity. AT 
is the temperature drop across the mantle minus the adi¬ 
abatic temperature change. We use the mantle viscosity 
calculated with (T m ) and (P), which is characteristic of 
the whole mantle layer. This choice dictates the behav¬ 
ior of our model, as will be discussed later. The viscosity 
parameterization is described in Section 2.4. 

Two other parameters are derived from the thermal 
model. The areal spreading ra te is t h e rate at w hich 
new ocean crust is being created. [McGovern fe Schubert I 
(|1989f ) parameterized the areal spreading rate using the 
current volume of the oceans and the present day heat 
flux, whic h are u nconstrained for exoplanets. Instead, we 
follow ISandu et al~l (| 2011 lh who relate the areal spread¬ 
ing rate to the convective velocity u c ) and the length of 
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Table 2 

Nominal Thermal Model Parameters. 


par am. 

units 

Upper Mantle 

Lower Mantle 

Core 

Rttcrit 


1100 

0.28-Ra 0 ' 21 

— 

a 

xlO -5 K” 1 

2 

1 

— 

k 

W m” 1 K- 1 

4.2 

4.2 

— 

K, 

m 2 s -1 

to -6 

10" 6 

— 

c p 

J kg" 1 K- 1 

1200 

1200 

840 

p 

kg m —3 

3300 

{pm) a 

8400 

a See Table 1 


Table 3 

Water Cycle Parameters. 

par am. 

name 

value 

Xr 

regassing efficiency 

0.03 

Xd 

degassing efficiency 

0.02 

fbas 

hydrated basalt fraction 

0.03 

Pbas 

basalt density (kg m —3 ) 

3000 

^ridge 

Mid-ocean ridge length 

1.5x(2ttHp) 

K 

Solidus depression constant (K wt% — 

43 

7 

Solidus depression coeffiient 

0.75 

6 

Melt fraction exponent 

1.5 

Dh 2 o 

Silicate/melt partition coefficient 

0.01 


the spreading centers ( L r id ge ) where ocean crust is cre¬ 
ated: 

S = ‘2L r j i( ig e , u JC (7) 

The convective velocity is determined by the convec- 

tive layer overturn ti me from boundary layer theory 
dSchubert, et al. ifeOOlL ch. 8 ) using the equation: 

_ 5.38 k{R p - R c ) to ^ 

u c — s-2 f°) 

Oy, 

We parameterize the length of the mid-ocean ridges as 
1.5 times the planetary circumference. This parameteri¬ 
zation is chosen to give the present day mid-ocean ridge 
(MOR) length on the Earth of ~60,000 km. We describe 
results using smaller L r id ge values in Section 4. 

Values for the parameters used in the thermal evolu¬ 
tion model are given in Table 2 . The value of R a rr it for 
the mantle is taken from lSchubert et al. I (120011) . We use 
constant values for the heat capacity, thermal conductiv¬ 
ity, thermal expansivity and thermal diffusivity, although 
these parameters are all known to be pressure-dependent. 


2.3. Volatile Evolution Model 

Vola tile evolution models harken back to 
iMcGovern fe Schubert I (|1989f ). repeated with vari¬ 
ations by many others. The volatile evolution model 
involves calculation of outgassing and ingassing rates for 
water based on mantle melting and surface hydration. 
The rate of change of the water abundance in the mantle 
is given by combining the ingassing and outgassing 
rates: 

dM H 2 o,m _ „ /rA 

j, — T'ingas Toutgas 


rameterized as: 


T'outgas — PrrijV^m^ (10) 

where p m , v is the density of volatiles in the mantle,, d m is 
the depth of melting, and S is the areal sprea ding rate of 
the mid-ocean ridges. In IMcGovern fc Schubert I dl989t) . 
d m is kept at a constant value of 100 km. The ingassing 
rate is parameterized as: 


T ingas 


fbas Pbasdbas $ X.r 


(ii) 


where fbas is the mass fraction of volatiles in the hy¬ 
drated basalt layer, pb as is the density of basalt, db as is 
the average thickness of the basalt (held constant at 5 
km) and \ r is an efficiency factor reflecting the incom¬ 
plete transport of the water in the hydrated basalt layer 
into the deep mantle. With all parameters except S and 
Pm,v held constant in equations (10) and (11), we found 
that all planets necessarily reached a steady state (i.e., 
= 0 )> where the mantle water abundance is 
given by setting r ingas equal to r outg as , and solving: 


_ ^^H 2 0,m _ fbasPbasdbasR\r 

V m d m S 


( 12 ) 


MH 2 0,m 


fbas PbasdbctsXrVm 


(13) 


While a case may be made that the Earth is in a steady 
state, there is little reason to suppose that this is a 
necessary condition for all exoplanets experiencing plate 
tectonics. W e therefore follow here the volatile evolu¬ 
tion model of lSandu et al. I (1201 H I. described briefly be¬ 
low. In this model, the depths of melting and hydration 
are not held constant, but vary based on local tempera¬ 
tures. We found that in these models, steady-state was 
rarely achieved. Parameters used in the volatile evolu¬ 
tion model are given in Table 3. The outgassing rate is 
given by: 

Poutgas = Pm{Rmelt) (-^-melt)Rmelt^Xd (M) 

where ( F me i t ) and (X me i f ) are the average fraction of 
melting and average abundance of water in the melt over 
the melt layer thickness, D me i t , S is the areal spreading 
rate and Xd is the degassing efficiency, which accounts for 
i ncomplete transport o f water to the surface. Whereas 
IMcGovern fe Schubert I d l989j) used a c onstant value for 
the melt layer thickness. iSandu et al. I (12011H used the 
mantle thermal profile and the peridotite solidus to deter¬ 
mine the melt layer thickness. The thermal profile used 
is composed of the conductive upper thermal boundary 
layer and the upper mantle adiabat, and the intersection 
of this profile with the hydrated solidus curve for peri¬ 
dotite determines where melt forms (see Fig. 1). Water 
dissolved in silicates lowers their solidus (the temperature 
at which partial melting begins), and the water partitions 
pr eferent ially into the melt. Using the parameterization 
of lKatz et al. I (I2003H for wet melting, the solidus depres¬ 
sion is given by: 

Rsol,wet — Tsoi,dry XT]j 2 o = T so l,dry R ^rnelt (1^) 


This equation is solved simultaneously with the heat 
t ransfer eq uation at eac h time. In the simplest form of 
IMcGovern fe Schubert I (|1989h . the outgassing rate is pa¬ 


where K and 7 are empirically determined constants for 
peridotite (see Table 3). The melt fraction and water 
abundance in the melt are determined where the mantle 
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thermal profile is above the wet solidus temperature and 
are given by: 


Emelt — 


T - 71 


n S 


sol,wet 


Tliq^dry Tsoi,dry 


(16) 


■Amelt — 


X H 2 0 ,m 


Dh 2 0 + Fmelt{ 1 — Dh 2 o) 


(17) 


where Xu 2 o,m is the water mass fraction in the mantle, 
Dh 2 o is the silicate/melt partition coefficient, and 6 is 
an empirically determined exponent. The dry liquidus 
and so lidus equations are taken from Zhang fc Herzbergl 


(11994T ) and iHirschmann et al~l f 2009h . respectively. The 
melt fraction and water fraction are averaged over the 
melt zone thickness at each timestep for use in equation 
(9). 

The return of water from the surface back into the 
mantle through ocean plate subduction is described by: 


Tingas — fbas Pbas Fh ydr S\ r 


(18) 


where fbas is the mass fraction of water in a hydrated 
basalt layer of thickness Dhydn Pbas is the density of 
basalt, and \r is the regassing efficiency, which accounts 
for imperfect return of water to the mantle. The hy¬ 
drated layer thickness is measured from the surface, down 
to the depth at which the temperature reaches the stabil¬ 
ity boundary of serpentinite (i.e., serpentine is not stable 
at lower depths). The upper thermal boundary layer has 
a linear conductive temperature profile (see eq. 5), where 
the temperature change with depth is ( T u —T s )/S u . The 
depth to the serpentine stability temperature is therefore 
given by: 


D 


hydr 


T 1 — T 1 T — T 

_ ^ serp i s _ serp s 


Tu~T s 


Qm 


(19) 


where T serp is the highest temperature of serpentine sta- 
bilit y at pressures be low ~ 3 GPa and is ~ 973 K 
ljUlmer fe Trommsdorffl 119951) . The hydrated layer is 
necessarily smaller than the thermal boundary layer, and 
is further restricted to hold no more water than is present 
at the surface in a given instant in order to maintain wa¬ 
ter mass-bal ance. We use the v alues for \r and \d de¬ 
termined bv lSandu et al. 1 (1201111 . 

A final word about parameterized volatile evolution 
models: As noted by a reviewer, the equations described 
above assume that water transported into the mantle is 
instantaneously transported throughout the mantle and 
available for outgassing. In the real world, of course, the 
mantle is not homogeneous and the spreading zone melt 
centers will likely become dehydrated before they can be 
replenished by advection from subducted water. There¬ 
fore outgassing rates calculated here are upper limits on 
the true values. 


2.4. Viscosity 

Water dissolved in silicate miner als su c h as o livin e re¬ 
duces the mineral’s strength (e.g., iChopra fc Paterson! 
EMU), and therefore its viscosity. In experimental liter¬ 
ature, the water dependence of the viscosity is parame¬ 
terized via the water fugacity / h 2 o , which depends on 
temperature and pressure. The rheological or constitu¬ 
tive law relating stress r and strain rate e for olivine is 


Table 4 

Viscosity Parameters. 


parameter 

Olivine (wet) 

Perovskite(dry) 

VO (Pa s) 

1.08 x 10 25 

1 x 10 21 

r 

1.0 


E a (kj mole -1 ) 

335 

300 

Va (cm 3 mole -1 ) 

4.0 

2.5 

Rgas (J mole -1 K -1 

) 

8.314 


then given by: 

e = A crp T n d- p lnr H20 exp (20) 

where n = 1 for diffusion creep, A crp is a material param¬ 
eter, d is the grain size in microns, E a is the activation 
energy, V a is the activation volume, and R gas is the ideal 
gas constant. This equation shows that the response to 
stress depends on the pressure and temperature, as well 
as the water fugacity ( Jh 2 o ) in a non-linear way. The 
viscosity is then derived from the constitutive law: 

r i e f f = ^ ( 21 ) 

We combine the constant parameters A crp and the grain 
size into a normalization factor etao to arrive at the ef¬ 
fective viscosity. The form of the effective viscosity is 
then: 


Veff = Vof H ^o ex P 


E a 1 

Rgas T 



1 ( PVa 

Rgas T 


PrefVg , 
T re f 


( 22 ) 


We normalize so that 7?(T re / = 1600, Xh 2 o= 500 ppm, 
P re f = 0)= 10 21 Pa s, which gives reasonable viscosities 
for the Earth. Param eters are given in Table 4 ; and are 
taken primarily from iHirth fe Kohlstedt I (1200311 for wet 
diffusion in olivine. These authors measured values of r 
of ~ 0.7 — 1.2 on polycrystalline olivine samples. Recent 
work on Si diffusion in olivine and single-crystal deforma¬ 
tion experiments suggests that the dependenc e on water 
abundance may be much l ower (r = 0 — 0.33) (Fci et al. I 
120131 [Girard et aT7~l I2013T) . However, more work needs 
to be done to reconcile these experiments with the poly¬ 
crystalline experiments. In this paper we use results from 
the earlier studies for the nominal models, but we also 
explore the effect of different r values on our findings. 

To convert from mass fraction water in t he ma ntle to 
fugacity / h 2 o , we use the formulation of lLi et al~l (l2008f ) . 
which is an empirical relationship between the water fu¬ 
gacity and the concentration of water in olivine ( Coh , 
atomic H/10 6 Si): 

l n fH 2 o = Co + c\lnCoH + C2ln 2 CoH + c^Iti^Coh (23) 

where Co = — 7.9859,Ci = 4.3559, C 2 = —0.5742, and 
C 3 = 0.0227. This relationship is derived from olivine 
solubility data between ~ 1373 — 1600 K and is only 
strictly valid within that temperature range. However, 
the relationship varies only marginally with temperature 
for a wide range of Coh and / h 2 o, so we apply it to 
the whole temperature range considered here for lack 
of other available data. It is straightforward to convert 
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Table 5 

Initial Temperature Parameters. 


Mass 

(M®) 

Cm 

(Trn,i) 

(K) 

T 

-*■ C,l 

(K) 

1 

1.19 

3000 

3810 

2 

1.32 

3330 

4630 

3 

1.44 

3630 

5350 

5 

1.64 

4130 

6640 


from the mantle water mass fraction to the concentration 
in Coh using the molecular weights of water and olivine. 

2.5. Initial parameters 

The initial temperature parameters for each of the 
super-Earth models are given in Table 5. Previ¬ 
ous parameterized convection models have shown that 
initial temperatures do not significantly affect the 
the rmal evolution beyond a few hundred Myr (se e 

e.g. lSchubert et ah~lll98d : lMcGovern fe Schubert 111 989h . 

Therefore, the present models are all started with the 
same initial mantle potential temperature T p , which is 
the temperature of the mantle adiabat extrapolated to 
the surface. The initial mantle potential temperature is 
set to 2520 K for all models, which is equivalent to an 
average mantle temperature of ~ 3000 K for the Earth. 
The average mantle temperatures are calculated accord¬ 
ingly for all super-Earth models. In models described in 
a later section, which include core evolution, we assume 
an initial temperature contrast across the CMB of 100 
K. Therefore, the initial core temperature is set equal 
to the temperature extrapolated along the mantle adia¬ 
bat to the CMB plus 100 K. We describe the effect of 
different initial temperatures in Section 4. 

3. RESULTS 

For our nominal model, we focus on the evolution of 
temperature and water abundances in both the mantle 
and on the surface. In the following sub-section, we in¬ 
troduce a second model which includes core-cooling and 
uses a different characteristic mantle viscosity for the up¬ 
per mantle that produces significantly different results. 

3.1. Nominal Model - Single Layer Convection 

Figures 2 and 3 show results for the nominal model 
using parameters from Table 3 and a mantle water 
abundance of ~1400 ppm (4 OM for 1 M®). Figure 2 
shows the evolution of the mantle potential temperature, 
viscosity and water abundance for the nominal model 
using two different viscosity parameterizations. For the 
curves shown in red, the activation volume is set to 0, so 
the viscosity is pressure-independent. The blue curves 
include a non-zero activation volume (see Table 4), so 
the viscosity depends on pressure as well as temperature 
and water fugacity. The viscosity is calculated with the 
average mantle temperature (T m ) and the mid-mantle 
pressure. 

For the pressure-independent viscosity, Fig. 2a seems 
to indicate that the large planets cool off more rapidly 
than the smaller planets. The 5 M® planet has a final 
potential temperature ~240 K lower than that of the 
1 M® planet. This is counter to standard intuition: 
large planets should cool off more slowly than small 





Figure 2. Nominal model with single layer convection, (a) man¬ 
tle potential temperature, (b) mantle viscosity, and (c) mantle wa¬ 
ter mass fraction. Red lines represent calculations done with a 
pressure-independent viscosity, blue lines represent the use of a 
pressure-dependent viscosity. Line styles indicate planet mass ac¬ 
cording to the legend. 
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Figure 3. Surface water abundance for the nominal models. The 
water abundance is equivalent to 4 ocean masses of water in the 
Earth-sized planet. 

planets due to their smaller surface-area/volume ratios. 
However, note that these are potential temperatures 
(i.e, the mantle temperature extrapolated to the surface 
along an adiabat). The average mantle temperature of 
the 5 planet is in fact hotter than the 1 M®, but by 
only ~80 K. The adiabatic gradient is much steeper for 
the large planets due to dependence on gravity, so the 
near-surface temperatures are therefore much lower. As 
seen in Fig. 2b, the hotter average mantle temperatures 
of the large planets results in lower mantle viscosities 
and therefore more rapid cooling. 

The lower near-surface temperatures of larger planets 
reduces the degree of melting, which strongly affects 
their water cycles, shown in Fig. 2c. There is a sharp 
decrease in mantle water abundance in the first 2 Gyr 
due to rapid early outgassing for all planets. As mantle 
temperature drops, near-surface temperatures drop 
below the solidus temperature, which halts melting and 
outgassing. This is followed by a rapid ingassing period. 
Ingassing is limited by the mass of the water at the 
surface and the thickness of the hydrated surface layer, 
which is regulated by the surface heat flux. All but the 
smallest planet have ingassed nearly all of their water 
by ~2 Gyr. The deep water cycles of these planets have 
effectively ceased. However, it should be noted that for 
all planets, this leaves a small residual surface reservoir 
of water that effectively cannot be lost to the mantle. 

For the pressure-dependent viscosity calculations (blue 
curves) both the average and potential temperatures of 
the smaller planets cool more quickly than the larger 
planets (see Fig. 2a). In fact, the largest planets 
initially heat up substantially before beginning to cool, 
so near-surface melting persists for much longer than in 
the pressure-independent case. The final temperature 
of the 5 M® planet is ~1000 K higher than that of the 
1 M® planet. The slow cooling of the planets can be 
attributed to the substantially larger viscosities in this 
model, due primarily to the pressure dependence of the 
viscosity (see Fig. 2b). The changes in water abundance 
shown in Fig. 2c are much more gradual than with 
pressure-independent viscosities. There is a rapid early 
outgassing phase only for the smallest planet, which 


has the lowest mantle viscosity, but this outgassing is 
much less complete than for the pressure-independent 
case. The larger planets show steady but very gradual 
outgassing, as their temperatures increase and their 
viscosities drop. In fact, the 5 M® planet has delayed 
onset of outgassing by ~1 Gyr, due to a very large 
thermal boundary and low surface heat flux, both of 
which can be attributed to the large initial mantle 
viscosity. 

The activation volume that we use here for the 
olivine viscosity is 4 cm 3 mol -1 , which is derived from 
experimental data. However, the activation volume 
for other sili cates has been shown to decrease with 
pressure (iStamenkovic et al. 1120111) . In their thermal 
models, IStamenkovic et al. I (|2012f) use activation vol¬ 
umes for pressures at the planet’s CMB (2.5 cm -3 
mol -1 for Earth). Therefore, a slightly lower activation 
volume may be more appropriate, particularly for 
the larger planets. Intermediate values of V a give 
mantle temperatures intermediate to those shown in 
Fig. 2a for the pressure-dependent viscosities. The 
mantle water abundance drops less precipitously in 
early times, and ingassing is much slower than for the 
pressure-independent case, but more complete than for 
the pressure-dependent model shown in Fig. 2b. The 
pressure-dependence of the viscosity is therefore an 
important parameter to include in the models. The 
value of activation volume will affect the final volumes 
of the surface oceans and the mantle temperature. 

Figure 3 shows the surface abundance of water for the 
nominal model for each of the super-Earths. This figure 
is the surface corollary to Fig. 2c. We show only the 
pressure-dependent models for clarity. Note that while 
the relative abundance of water is the same for each 
planet, the total planetary water mass is 4 (8, 12, 20) 
ocean masses for the 1 (2, 3, 5 M®) planet. The time- 
dependent behaviors of the planets do not scale simply 
with planet mass. The 1 M® planet has a much more 
significant outgassing phase at early times, followed by 
ingassing. The larger planets show much more gradual 
outgassing, with delayed onset of outgassing for the 5 
M® planet. Due to the delayed outgassing, the 5 M® 
planet has less surface water than the 3 M® planet 
for most of their lifetimes. However, the 3 and 5 M® 
planets have roughly equal surface water abundances at 
10 Gyr. For smaller values of the activation volume, the 
surface oceans will be substantially larger. 


3.2. Boundary Layer Convection 

Volatile evolution models have typically assumed 
single-layer convection, which has been shown to work 
well for the Earth. However, thermal evolution models 
that neglect volatile evolution often assume that mantle 
convection is controlled by boundary layer instability, in 
which convective instability is determined by the local 
conditions at the boundary layers rather than the man¬ 
tle as a whole. We will now describe how these mod¬ 
els differ from the nominal model described above. For 
the boundary layer models, we assume mixed heating, 
with a non-zero heat flux from the core (taken here to be 
both solid and isothermal) and a conductive lower bound¬ 
ary layer. Here we are less interested in reproducing the 
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Earth than in exploring the behavior of the models for 
different planet masses. When using the same physically- 
plausible parameters here, we show that the two types of 
models give fundamentally different results, due to the 
choice of characteristic mantle viscosity. 

For the boundary layer models, a second heat transfer 
equation is needed for the core: 


little variability in the results for surface oceans when 
holding a constant. 

Results for the boundary layer convection model are 
discussed in the following subsection. Afterwards, we 
discuss how parameters affect the two different models. 

3.2.1. Results - Boundary Layer Convection Model 


dT 

PcCp cYc 7~ = A c q c 
at 


(24) 


where variables are as in eq. ( 1 ), but defined for the core 
rather than the mantle. The core is isothermal, and so 
is characterized by a single temperature T c . We neglect 
radioactive decay in the core, as well as latent heat due 
to possible core freeze-out. The heat flux out of the core 
is given by: 


Qc = 


k 


0 T c - T,) 
Sc 


(25) 


where S c is the lower thermal boundary layer, T c is the 
core temperature, and Tj is the mantle temperature at 
the top of the lower thermal boundary layer (see Fig. 1). 
We use a local Rayleigh number to define the thickness of 
the lower boundary layer. The boundary layer thickness 
is determined by setting it equal to its critical thickness, 
the point at which it becomes unstable to convection. 


1 

/ Kr](T c , P cr nb)kla cr it,l \ 3 
v gap m (T c -Ti ) ) 


where 77 (T c , P cm b) is the viscosity of the lower bound¬ 
ary layer. For the lower boundary layer, we use a per- 
ovskite rheology. The choice of rheologies will be dis¬ 
cussed more below. The critical Rayleigh number of the 
l ower bo unda ry layer is giv en by Ra cr i tj i = 0.28i?a°' 21 
(iDeschamps fc Sotin I |2000h . where Ra is the global 
Rayleigh number used in the definition of S u . The global 
Rayleigh number is defined here as in eq. ( 6 ), except for 
the characteristic viscosity. For the characteristic vis¬ 
cosity, we use the value defined by the temperature and 
pressure at the base of the upper thermal boundary layer. 
This has a signifcant effect on the outcome of the models 
as shown below. 

The viscosity of olivine is used for the upper thermal 
boundary layer, but olivine is not a stable phase at the 
pressures and temperatures of the lower thermal bound¬ 
aries. We theref ore use th e vis cosity law for perovskite 
derived bv IStamenkovic et al. 1 (|20111) for the lower ther¬ 
mal boundary. Parameter values are given in Table 4. 
The lower mantle of super-Earths likely consists of per¬ 
ovskite t ransitioning to post-perovskite for larger pl anets 
(see e.g. IValencia et al~1 120061 [Ta.cklev et al. 1120 lit) . No 
experimental data exists on the water-dependence of the 
viscosi ty of perovs ki te, so the value of r is set to zero. 
IStamenkovic et al~l (l 201 lh also give the dependence of 
the activation volume as a function of pressure. In the 
nominal models, we use a value of 2.5 cm 3 mol -1 , which 
is the value at the Earth’s CMB. 

We use a slightly smaller value of a for the lower 
therm al boundary layer (see Table 2), which is taken 
from iTacklev et al. I (1201311 for the Earth’s CMB. 
IStamenkovic fe Breuer I ()2014l ) found that scaling a with 
planetary mass was as important as scaling p m to the 
calculation of planetary temperatures. However, we find 


Results are shown in Figure 4. Core temperatures are 
normalized to their initial values, because they vary by 
~3000 K. The upper boundary layer is calculated using 
the pressure-dependent olivine viscosity using temper¬ 
ature and pressure at the base of the boundary layer. 
The lower boundary is calculated using t he pressur e- 
depen dent perovskite viscosity from IStamenkovic et al.~l 
( I 20 TII) using temperatures and pressures at the CMB. 
This viscosity does not depend on the mantle water 
abundance. 

Mantle potential temperatures decrease more rapidly 
and significantly here than for the nominal model. All 
planets have nearly identical potential temperature 
evolution, in contrast to the pressure-dependent results 
shown in Fig. 2a. The upper thermal boundary layer 
thickness is comparable for all planets, so the upper 
mantle viscosity used here is only weakly dependent on 
pressure (blue curves, Fig. 4c). In contrast, the viscosity 
of the lower thermal boundary layer is strongly depen¬ 
dent on pressure, and therefore the core temperature 
(Fig. 4b) is highly dependent on planet mass. The 1 
M® planet’s core cools by ~35%, whereas the 5 M® 
planet’s core cools only ~5%. 

Figure 4d shows the surface water inventories for 
the boundary layer model, with 4 ( 8 , 12, 20) OM of 
initial water. In comparison to Figure 3, it is obvious 
that significantly more water is outgassed from the 
mantle here. The 5 M® planet has a peak surface water 
abundance of ~8 OM, in comparison to ~1.6 in Fig. 3. 
However, the residence time is much shorter. All of the 
planet’s lose most of their surface water inventory back 
into the mantle by ~4.5 Gyr. After complete ingassing, 
the 1 M® planet has the largest remaining surface water 
abundance, with ~0.2 OM at 10 Gyr. All planets have 
significantly lower surface water abundances at 10 Gyr 
for the boundary layer models than for the nominal 
models shown in Figure 3. For the boundary layer 
models, we find little difference between the pressure- 
dependent and pressure-independent viscosities, except 
for the cooling of the core. The peak surface water 
abundances remain essentially unchanged, but the pres¬ 
ence of surface water does persist for ~0.5 Gyr longer 
in the pressure-independent case. This is not surprising, 
since the pressures of the upper thermal boundary layer 
are very low and do not signficantly affect the viscosities. 


3.3. Dependence on parameters 

Many of the parameters used in the thermal and 
volatile evolution models are poorly constrained for the 
Earth, much less super-Earths. In the following section, 
we exam the effect of varying several of these parame¬ 
ters on the results of both the nominal and the bound¬ 
ary la yer convection models. We refer the reader to 
IStamenkovic fe Breuer I (|2014f) for an analysis of the ef¬ 
fect of /?, a m , K m , and p m on thermal evolution models. 

































Oceans on Super-Earths 


9 





Figure 4. Results for nominal models with boundary layer convection heating, (a) mantle potential temperature, (b) core temperature 
normalized to initial core temperature, (c) thermal boundary layer viscosities (upper = blue, lower = red), (d) surface water abundance in 
ocean masses. Line styles indicate planet mass according to the legend. 


3.3.1. Fugacity coefficient 

Figure 5 explores the effect of the fugacity coefficient r 
on the results of the previous models for planets of 1 and 
5 M©. Other planet masses are not shown for clarity. 
This fig ure shows r val u es of 0.7, 1.0 (nominal), and 
1.2. iHirth &; KohlstedTI (12003H give values of 0.7 — 1.0 
for wet diffusion creep of olivine and r = 1.2 for wet 
dislocation creep of olivine, which is why we chose these 
values. It should be noted that the viscosities were not 
renormalized with the change of r values. Therefore the 
variances reflect the absolute changes in viscosity. 

For the single layer convection model (Fig. 5a), the 
initial outgassing phase is significantly stronger for 
r=1.2 for the 1 M© planet, but the final abundance 
is nearly the same as the nominal value. In contrast, 
outgassing is slightly delayed for r = 0.7 but the final 
water abundance is ~0.2 OM larger. For the 5 M© 
planet, the larger r value increases the amount of water 
outgassed, whereas the lower r value both further 
delays outgassing and limits the total amount of water 


outgassed significantly. 

For the boundary layer convection model (Fig. 5b), 
the higher r value increases the amount of outgassing 
for both the 1 and 5 M© planets and shifts the peak of 
outgassing to slightly earlier times. However, ingassing 
also occurs more rapidly, so the oceans persist only 
until ~4 Gyr. The lower r value reduces the amount 
of outgassing for both planets, and causes the surface 
water to persist for longer. For the 1 M© planet and r 
= 0.7, the water abundance is relatively constant over 
the planet’s lifetime. There is about 0.4 OM of water 
remaining on the surface after 10 Gyr. 

3.3.2. Total water abundance 

Figure 6 shows the surface water abundances using 
different initial mantle water abundances. The surface 
water abundances appear to be a fairly straightforward 
function of water abundance. Models with larger 
water abundance have larger surface inventories. The 
outgassing occurs earlier for the single layer convection 
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Figure 5. Surface water abundances for (a) the nominal single layer convection model and (b) boundary layer convection models using 
different values of the fugacity exponent r. Blue lines are for 1 M®, red lines for 5 M®. 




Figure 6. Surface water abundances for (a) the nominal single layer convection model and (b) boundary layer convection models for 
different initial water abundances. Total water abundances are equivalent to 2, 4 and 6 ocean masses of water for the Earth-mass planet. 
Abundances for the larger planets are the same in terms of mantle mass fraction of water. Colors indicate planet mass (blue 1 M®, red 
5 M 0 ), and line styles indicate water abundance. 


models (Fig. 6a), and oceans persist later for the 
boundary layer convection models (Fig. 6b). The 
effect of varying initial water abundances produces 
results similar to changes in the fugacity coefficient 
shown in Fig. 5. However, one major difference to 
note is that although both larger water abundances 
and larger fugacity coefficients increase outgassing, 
the persistence of the oceans differs. Surface oceans 
persist longer (til ~8 Gyr for the 5 M® planet) for en¬ 
hanced water, whereas increasing the fugacity coefficient 
shortened the ocean lifetime. Another thing to note 
is that the surface water abundances for the boundary 
layer model with 2 (10) OM of water are nearly identical. 


3.3.3. Initial mantle temperature 


Although models for the Eart h show limited sen sitiv¬ 
ity to in itial temperature ( e.g. iMcGovern fc Schubertl 
(jl989f) . iTaiika fc Matsui I (1 199T b this appears not 
to be the case for the super-Earth models. Figure 
7 compares results for the single layer model for the 
nominal starting mantle temperature of 2520 K, 2000 
K, and 3000 K for the 1 (blue) and 5 (red) M® planets. 
Mantle temperatures for the smaller planets converge 
within ~ 4 Gyr on the nominal results (Fig 2a) but 
for the hotter starting temperature the 5 M® planet 
remains persistently hotter throughout its lifetime. The 
hotter initial temperatures effect the water cycle for all 
of the planets. Although the 1 M® planet converges to 
nearly the same temperatures, the initially hotter planet 
outgasses 3x more water within the first 500 Myr. The 
water is gradually ingassed over the planet’s lifetime, 
but the final abundance of water remains slightly larger 
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Figure 7. Surface water for the single layer convection model. 
Solid lines show the nominal model (also shown in Fig. 2 and 
3), compared with models with either a higher (3000 K) or lower 
(2000 K) starting mantle potential temperature. Line styles indi¬ 
cate temperature according to the legend. Colors indicate planet 
mass (blue 1 M®, red 5M®). 

than for the colder starting planet. The 5 M® planet 
begins outgassing much more rapidly, and continues 
steadily outgassing for its lifetime. It ends with ~2 OM 
of water on the surface. Lower initial temperatures, 
which are more widely used in the literatur e (e.g. 
iStamenkovic et al~ll2012t iNoack fe Breuer 1120l.'lll . delay 
and reduced the degree of outgassing, particularly for 
larger planets. For an initial potential temperature of 
2000 K, the 5 M® planet does not begin degassing until 
~4 Gyr. Given th at su r face water is likely necessary for 
plate tectonics (jKorenaga I l2010bi ). these planets may 
not experience plate tectonics at all, until very late in 
their lifetimes. These planets will likely start in a stag¬ 
nant lid mode, which we have not attempted to model 
here. However, it is likely that the stagnant lid would 
allow the planet to heat up earlier so that formation of 
the oceans may not be as delayed as shown here. The 
boundary layer models show very limited dependence on 
the initial temperature and so are not shown here. For 
an initial starting temperature of 3000 K, the planets 
evolve at virtually the same temperatures, and the 
surface water abundances are only slightly enhanced. 
The different initial water abundances shown in Fig. 6 
had a far larger impact on that model’s results. 


3.3.4. Mid-ocean ridge length 

While we use a parameterization that produces 
the Earth’s present day mid-ocean ridge length, we 
note that the length of the ridges on Earth may have 
changed throughout time. However, results for both 
models change only slighlty with different values for 
Lridge ■ For the single layer model, the temperatures 
with variable L r id ge do not change, and surface water 
abundances only slightly decrease. The effect on the 
water abundance for the boundary layer models is 
moderate, as shown in Figure 8. We show results 
for models using lower values of L r id ge , which would 
correspond with slower plate growth. The effect on 



time (Myrs) 


Figure 8. Surface water abundances for the boundary layer con¬ 
vection model. Solid lines show the nominal model (also shown 
in Fig. 4), compared with models using either L r id ge = 100% 
(dash) or 50% (dash-dot), respectively of the planetary circum¬ 
ference. Colors indicate planet mass (blue 1 M®, red 5M®). The 
nominal model uses L ri ^ ge = 150% of the planetary circumference. 

temperature is minor for both the 1 and 5 M® planet, 
so we do not show it here. For the surface water 
abundances, the lower L r id ge values reduce the peak 
surface water abundance and significantly prolong the 
ingassing of the surface water after the initial outgassing 
phase. The 5 M® planet has ~0.5 OM of surface water 
remaining at 10 Gyr, whereas the 1 M® planet has 
~0.8 OM of surface water remaining for an L r id ge = 50%. 


3.3.5. Other parameters 

Another parameter that can significantly affect the 
model results is the viscosity parameterization. We 
have chosen here to normalize the viscosity to 10 22 Pa 
s for a reference state of 1600 K, 0 Pa, and 500 ppm 
water. This gives a reasonable viscosity for the Earth’s 
mantle. Howe ver, t here are other choices that could 
be made. iSandu et al. I (l201ll) tuned their model to 
match the present day Earth’s viscosity and heat flux, 
by normalizing their viscosity to 2.3 x 10 21 Pa s at 2300 
K and 500 ppm of water. Using this reference value 
for the mid-mantle pressure, we get 770 = 2.5 x 10 25 
Pa s, roughly an order of magnitude lower than the 
value used here. For the single layer model, using this 
normalization factor results in mantle temperatures 
cooler by ~100-150 K, and surface water abundances 
that resemble the hotter starting temperature in Fig. 7. 

The abundance of radioactive elements affects planets 
in both models. A recent paper shows that older planets, 
those formed early in the galaxy’s lifetime, will have 
much lower abundances of radioactive elements, whereas 
younger planets may h ave u p to 7 time s more radioactive 
heat production (jFrank et al. I 12014 1. For the single 
layer models, increasing the uranium abundance by 
a factor of two increases the mantle temperatures by 
100-200 K. In particular, all but the 1 M® planet heat 
up within the first 2 Gyr, rather than cooling. The 5 
M® planet reaches a peak temperature of ~3200 K, 
compared to ~3000 K for the nominal results. The 
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surface water abundances are increased slightly, but to a 
lesser extent than seen for a raise of initial temperature 
(Fig. 7). The mantle temperatures of the models heated 
from below increase by ~100 K. The peak surface water 
abundances do not change from the nominal results, 
but the ingassing of the water back to the mantle takes 
a longer amount of time. The 1 Af® has substantial 
surface water until ~ 4 Gyr, whereas the 5 M e planet’s 
surface water persists to ~ 6 Gyr, compared to ~4.5 
Gyr for the nominal models. 

4. DISCUSSION OF RESULTS 

4.1. Stagnant lid regime 

Many terrestrial planets, such as Mars and Venus, do 
not experience plate tectonics, but are in a stagnant- 
lid regime. For these planets, a thick lithosphere insu¬ 
lates the convecting portion of the mantle from the sur¬ 
face. Stagnant lids develop when planets are too cool to 
maintain convection across the entire silicate layer. It 
is likely that most planets can transition between plate 
tecto nics and the stagnant or sluggish lid regimes (iSleeo I 
12000( 1. However, the mechanisms by which this transi¬ 
tion occurs are poorly understood. We can speculate, 
however, on how the transition would effect the models 
described. The formation of a stagnant lid insulates the 
mantle and allows the temperature to increase, which 
would lower the viscosity and increase convective vigor. 
Regassing would halt in this regime, as there is no trans¬ 
port mech anism f or water to return to the mantle (see 
e.g. iMorschhauser et al. 112011 1. For all of the planets in 
the boundary layer model, the formation of a stagnant 
lid would likely extend the lifetime of the surface oceans, 
by halting regassing. The planets could also heat up 
enough to re-initiate mantle melting and degassing. Fu¬ 
ture work will look at the effect of such a transition on 
the persistence of surface oceans. 

4.2. Steady state versus thermal evolution 

As discussed in section 2.3, us ing th e vo latile evolu¬ 
tion parameterization of lMcGovern fc Schubert! (Il989ll . 
planets were found to achieve a volatile steady state. 
For a planet with a large global water abundance, most 
of the water will be found at the surface of the planet 
for the planet’s lifetime. For planets with global water 
abundances smaller than the steady state value, all 
water will be trapped within the mantle. Note that 
I Cowan fe Abbot] (120141 1 considered a steady-state solu¬ 
tion to determine the water mass fraction necessary for 
a planet to be considered a water planet (>90% surface 
ocean coverage). However, in the models considered 
here, steady state was never reached. In fact, the results 
discussed above show widely divergent outcomes for the 
same planet based upon the characteristic upper mantle 
viscosity chosen. Once the upper mantle temperature 
drops below the peridotite solidus, the melt layer disap¬ 
pears and the volatile cycle effectively ceases, although 
slow ingassing may continue until the surface becomes 
depleted of all water. 


4.3. Ocean depths on Super-Earths 


We have focused in the discussion of results on 
the surface water abundances of the different models 
because this is a potentially observable parameter, and 
one that has implications for the possible habitability 
of super-Earths. Concerns have been raised about the 
habitability of, or more importantly the ability of life 
to begin on, planets with global ocean coverage. Some 
fraction of continental surface, which provides a shallow 
water environment, may be necessary for life to begi n 
and fo r the evolution of complex life. iCowan fc Abbot! 
(12014( 1 derive a maximum ocean depth that separates 
planets with continents from totally water-covered 
planets using a crustal buoyancy model. They find that 
maximum ocean depth scales with surface gravity as 
dmax ~ 11.4 {g/g@) 1 km. Applying this to the results 
for the models above gives us the minimum surface 
area coverage on each planet. We show these results 
for both the single layer and boundary layer models 
in Figure 9. The single layer convection planets have 
slightly less than 50% areal coverage by oceans. The 5 
Mg planet has minimal ocean areal fraction until ~2 
Gyr, which would suggest that life would be difficult 
to begin on such a planet in its infancy. For boundary 
layer convection, the 3 and 5 Mg planets have areal 
coverages greater than 1, which indicates that they will 
have global oceans. These planets are not likely to have 
exposed continents. However, the 1 and 2 M® planets 
have significant continental area for both convection 
modes, which indicates that these planets may be best 
places to search for life. 


4.4. Additional planetary processes 

The models presented here are not comprehensive pa- 
rameterizations of processes that can affect the water 
budget of a planet’s surface. One particularly impor¬ 
tant proce ss is the loss of water f rom a planet’s at¬ 
mosphere ((Wordsworth fe Pierreliumbert 1 12013(1 . Our 
model is agnostic to the form that water takes at the 
surface (i.e., water, ice, steam, etc.), other than by 
constraining the surface temperature. However, the 
surface temperature will vary with stellar type, or¬ 
bital dista nce, atmospheric composition, and age. As 
IWordsworth &; Pierrehumbert I (I2013D show, loss of wa¬ 
ter vapor depends not only on the surface tempera¬ 
ture, but also on the CO 2 abundance of the atmo¬ 
sphere. Significant loss of water vapor on hot C02-rich 
planets would reduce the amount of regassing, so the 
mantle would become more dehydrated over time. We 
also negl e ct con t inen tal crust formation and weathering 
((Honing et al. I [2014T) . Continental crust effectively acts 
as an insulating barrier, but also serves as a sink for 
radioactive elements which are extracted from the man¬ 
tle. Weathering of continents and sedimentation rates, 
which are enhanced by living organisms, ca n affec t the 
regassing rate into the oceans. (Honing et al. I (12014(1 sug¬ 
gest that planets without life will evolve to have more 
surface water (therefore lower continental coverage), and 
dryer mantles than planets with life. Weathering also 
acts as climate control by stabilizing atmospheric CO 2 , 
which affects the re tention of wat er in the atmosphere 
as described above ((Abbot et al. ll20T2l f. Much further 
work needs to be done to truly understand the feedbacks 
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Figure 9. Minimum surface area coverage for the nominal models with a) single layer convection and b) b ound ary l ayer convec tion. The 
minimum surface area coverage is determined assuming the maximum ocean depth based on the scaling of ICowan He Abbot! (1201411 . Planets 
with greater than 90% ocean coverage are considered water planets. 


that will contribute to the persistence of oceans on super- 
Earths. 

5. SUMMARY 

We explored two different scaling parameterizations 
for plate tectonics planets using either single layer 
convection or boundary layer convection. Mantle tem¬ 
peratures are significantly hotter for the first model, and 
the surface water abundance lower, but more persistent. 
For many different parameters, smaller planets will 
have initially larger surface oceans, but lose them more 
rapidly than larger planets. Larger planets show delayed 
outgassing, which may compromise them as locations on 
which life can originate. The boundary layer convection 
model cools very rapidly due to vigorous convection 
driven by both an upper and a lower thermal boundary 
layer. The surface water abundances on the massive 
planets in the boundary layer model are extremely 
high and suggest that these planets will have limited, 
if any, continental coverage. However, these oceans 
persist for less than half of the planet’s lifetime. Upon 
complete ingassing of the oceans, the massive planets 
are effectively tectonically dead, and therefore unlikely 
to be habitable. 

Observations of rocky exoplanets in the 1-5 
range are already producing very accurate planet 
radii and mass determinations dDressing et al. II2014H . 
Many of these exoplanets have ages, determined from 
asteroseismic ages of their host stars. In the era of 
JWST and large ground-based telescopes, some of the 
nearest exoplanets will have atmospheric spectroscopy 
capable of distinguishing different geochemical regimes, 
like some of the extremes described here. Understanding 
the general features of the deep water cycle across 
rocky planets of different mass and structure might give 
us unique windows into their interior through its sin¬ 
gular effect on atmospheric and surface water abundance. 


We thank Li Zeng for helpful discussions and an anony¬ 
mous referee for a detailed review that greatly enhanced 


the quality of this paper. 
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